﻿using System;
using System.Collections.Generic;
using System.Diagnostics;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using Mapack;

namespace QProg
{
	public static class QProgUsual
	{
		#region FIelds

		public static int CountOfFuncValue { get; set; }
		public static int CountOfGradValue { get; set; }
		public static int CountOfGesseValue { get; set; }

		public static int Size { get; set; }
		public static double[] Dmin { get; set; }
		public static double[] Dmax { get; set; }
		public static double Deltaj0 { get; set; }
		public static double Alpha { get; set; }
		public static int J0 { get; set; }
		public static double TetaAlpha { get; set; }
		public static double Eps { get; set; }
		public static double TetaJ0 { get; set; }
		public static double Teta0 { get; set; }
		public static double Lj0 { get; set; }
		public static bool[] Optimal { get; set; }
		public static int I0 { get; set; }
		public static Matrix Poc { get; set; }
		public static Matrix Loc { get; set; }
		public static Matrix Teta { get; set; }
		public static Matrix D { get; set; }
		public static Matrix X { get; set; }
		public static Matrix C { get; set; }
		public static Matrix L { get; set; }
		public static double[] Delta { get; set; }
		public static Matrix Goc { get; set; }
		public static List<int> Joc { get; set; }
		public static Matrix A { get; set; }
		public static Matrix B { get; set; }
		#endregion

		public static double[] Grad(Matrix x)
		{
			double[] grad = new double[x.Rows];
			double t = 0;
			double expSum = 0;
			var temp = 0.0;
			for (int i = 0; i < grad.Length; i++)
			{
				t = (D.Submatrix(i, i, 0, D.Columns - 1) * x)[0, 0] + C[i, 0];

				for (int j = 0; j < A.Rows; j++)
				{
					expSum = 0.0;
					for (int k = 0; k < A.Columns; k++)
					{
						expSum += A[j, k] * X[k, 0];
					}
					temp += Math.Exp(expSum) * A[j, i];
				}
				grad[i] = t + temp;
				temp = 0;
			}
			return grad;
		}

		public static double[] ApproxGrad(Matrix x)
		{
			var delta = 1e-9;
			var e = new Matrix(x.Rows, 1);
			var grad = new double[x.Rows];
			for (int i = 0; i < x.Rows; i++)
			{
				for (int j = 0; j < x.Rows; j++)
				{
					e[j, 0] = 0;
				}
				e[i, 0] = 1;

				grad[i] = (CalculateF(x + e * delta) - CalculateF(x)) / delta;
			}
			return grad;
		}

		public static double Gesse(Matrix x, int k, int l)
		{
			var t = D[k, l];
			double temp = 0;
			double expSum = 0;
			for (int i = 0; i < A.Rows; i++)
			{
				expSum = 0.0;
				for (int j = 0; j < A.Columns; j++)
				{
					expSum += A[i, j] * X[j, 0];
				}
				temp += Math.Exp(expSum) * A[i, k] * A[i, l];
			}
			return t + temp;
		}

		public static double CalculateF()
		{
			var t = 0.5 * (X.Transpose() * D * X)[0, 0] + (C.Transpose() * X)[0, 0];
			double expSum = 0;
			for (int i = 0; i < A.Rows; i++)
			{
				expSum = 0;
				for (int j = 0; j < A.Columns; j++)
				{
					expSum += A[i, j] * X[j, 0];
				}
				t += Math.Exp(expSum);
			}
			return t;
		}
		public static double CalculateF(Matrix x)
		{
			var t = 0.5 * (x.Transpose() * D * x)[0, 0] + (C.Transpose() * x)[0, 0];
			double expSum = 0;
			for (int i = 0; i < A.Rows; i++)
			{
				expSum = 0;
				for (int j = 0; j < A.Columns; j++)
				{
					expSum += A[i, j] * x[j, 0];
				}
				t += Math.Exp(expSum);
			}
			return t;
		}



		public static void Solve(Matrix d, Matrix c, Matrix a, double[] dmin, double[] dmax, int size)
		{
			Stopwatch stopWatch = new Stopwatch();
			stopWatch.Start();
			D = d;
			C = c;
			A = a;
			Dmin = dmin;
			Dmax = dmax;
			Size = size;


			bool[] optimal = new bool[Size];
			var eta = 1.0;

			Joc = new List<int>();
			int countOfIterations = 0;
			L = new Matrix(Size, 1);
			X = new Matrix(Size, 1);
			Teta = new Matrix(Size, 1);

			GeneratorForMathematica.Generate(D, X, C, A, dmin, dmax, Size);

			Delta = Grad(X);




			var oldFuncValue = CalculateF();
			var newFuncValue = 0.0;

			for (int i = 0; i < Size; i++)
			{
				X[i, 0] = 0;
			}

			while (true)
			{
				CountOfGradValue++;
				Delta = Grad(X);
				var ax = ApproxGrad(X);

				for (int i = 0; i < Size; i++)
				{
					if (Math.Abs(Delta[i]) < 1e-8)
					{
						Delta[i] = 0;
					}
				}

				CountOfFuncValue++;
				newFuncValue = CalculateF();
				if (oldFuncValue < newFuncValue)
				{
					Console.WriteLine("Old: {0}", oldFuncValue);
					Console.WriteLine("New: {0}", newFuncValue);
					oldFuncValue = newFuncValue;
				}

				#region JocCheckness
				//if (Joc.Count > 0)
				//{
				//	for (var item = 0; item < Joc.Count; item++)
				//	{
				//		if (Math.Abs(Delta[Joc.ElementAt(item)]) > 1e-8)
				//		{
				//			//	if (Joc.Count > 0)
				//			//	{
				//			var index = item;
				//			//Запоминаем оставшиеся индексы
				//			var indexes = new List<int>();
				//			foreach (var it in Joc)
				//			{
				//				if (Joc.IndexOf(it) != index)
				//				{
				//					indexes.Add(Joc.IndexOf(it));
				//				}
				//			}
				//			//Console.WriteLine("Joc count: {0}", Joc.Count);
				//			var GocTemp = new Matrix(Goc.Rows - 1, Goc.Columns - 1);
				//			Joc.Remove(Joc.ElementAt(item));

				//			//Формируем новую матрицу Goc с выведенным индексом
				//			for (int i = 0; i < Goc.Rows - 1; i++)
				//				for (int j = 0; j < Goc.Columns - 1; j++)
				//				{
				//					GocTemp[i, j] = Goc[indexes[i], indexes[j]];
				//				}

				//			Matrix GocHorizontal = new Matrix(1, Joc.Count);
				//			for (int i = 0; i < Joc.Count; i++)
				//			{
				//				GocHorizontal[0, i] = Goc[index, indexes[i]];
				//			}

				//			Matrix GocVertical = new Matrix(Joc.Count, 1);
				//			for (int i = 0; i < Joc.Count; i++)
				//			{
				//				GocVertical[i, 0] = Goc[indexes[i], index];
				//			}

				//			Goc = GocTemp - GocVertical * GocHorizontal * (1 / Goc[index, index]);
				//			//	}
				//			//else
				//			//{
				//			//	var index = Joc.ElementAt(item);
				//			//	Joc.Remove(index);
				//			//	Goc = new Matrix(1, 1);
				//			//}
				//		}
				//	}
				//}
				#endregion


				countOfIterations++;
				if ((countOfIterations % Size) == 0)
				{
					Joc.Clear();
				}

				for (int i = 0; i < Size; i++)
				{
					if (
						!(
						((Math.Abs(X[i, 0] - dmin[i]) == 0) && Delta[i] >= 0)
						|| ((Math.Abs(X[i, 0] - dmax[i]) == 0) && Delta[i] <= 0)
						|| ((X[i, 0] > dmin[i] && X[i, 0] < dmax[i]) && Math.Abs(Delta[i]) < 1e-5))
						)
					{
						optimal[i] = false;
					}
					else
					{
						optimal[i] = true;
					}
				}
				if (optimal.All(z => z == true) || (newFuncValue > oldFuncValue + 1e-6))
				{

					stopWatch.Stop();
					TimeSpan ts = stopWatch.Elapsed;
					string elapsedTime = String.Format("{0:00}:{1:00}:{2:00}.{3:000}", ts.Hours, ts.Minutes, ts.Seconds, ts.Milliseconds);
					Console.WriteLine();
					Console.WriteLine("Классический");
					Console.WriteLine("RunTime: " + elapsedTime);
					Console.WriteLine("Количество вычислений значений функции: {0}",CountOfFuncValue);
					Console.WriteLine("Количество вычислений градиента функции: {0}", CountOfGradValue);
					Console.WriteLine("Количество вычислений гессиана функции: {0}", CountOfGesseValue);


					Console.WriteLine("Решение найдено");
					Console.WriteLine();
					//Console.WriteLine();
					//Console.WriteLine("Градиент");
					//var grad = Grad(X);
					//for (var i = 0; i < grad.Length; i++)
					//{
					//	Console.WriteLine("grad[{0}]={1}", i, grad[i]);
					//}
					Console.WriteLine("F[x] = {0}", CalculateF());
					Console.WriteLine("Количество итераций: {0}", countOfIterations);
					//for (int i = 0; i < Size; i++)
					//{
					//	Console.WriteLine("X[{0}]={1}", i, X[i, 0]);
					//}
					Console.ReadKey();
					break;
				}

				//Ищем Delta j0
				int k = 0;
				var flag = false;
				while ((optimal[k] != false) || (Joc.Contains(k)))
				{
					if (k == optimal.Length - 1)
					{
						Joc.Clear();
						flag = true;
						break;
					}
					k++;
				}
				if (flag)
				{
					k = 0;
					while (optimal[k] != false)
					{
						k++;
					}
				}
				flag = false;


				Deltaj0 = Delta[k];
				J0 = k;
				for (int i = 1; i < Size; i++)
				{
					if ((Math.Abs(Delta[i]) > Math.Abs(Deltaj0)) && (optimal[i] == false) && (Joc.IndexOf(i) == -1))
					{
						J0 = i;
						Deltaj0 = Delta[i];
					}
				}

				if ((J0 == 0) && (Joc.IndexOf(0) != -1))
				{
					Joc.Clear();
					k = 0;
					while (optimal[k] != false)
					{
						k++;
					}

					Deltaj0 = Delta[k];
					J0 = k;

					for (int i = 1; i < Size; i++)
					{
						if (Math.Abs(Delta[i]) > Math.Abs(Deltaj0) && optimal[i] == false)
						{
							J0 = i;
							Deltaj0 = Delta[i];
						}
					}

				}

				Lj0 = Math.Sign(Deltaj0);


				//Проверка условия равенства множества Joc пустому множеству
				if (Joc.Count == 0)
				{
					CountOfGesseValue++;
					Alpha = Gesse(X, J0, J0);

					if (Alpha < 0)
					{
						throw new Exception();
					}

					for (int i = 0; i < Size; i++)
					{
						L[i, 0] = 0;
					}
					L[J0, 0] = -Lj0;
					TetaAlpha = Math.Abs(Deltaj0) / Alpha;
					CountOfFuncValue+=2;
					eta = (CalculateF(X) - CalculateF(X + L * TetaAlpha)) / 2;

					//Func<double, double> f = (teta) => CalculateF(X + L*teta) ;
					//var min = GoldCut(f, -2, 2, 1e-15);
					//TetaAlpha = min;

					// Находим минимальное Teta
					Teta0 = TetaAlpha;


					if (Joc.Contains(J0))
					{
						throw new Exception();
					}
					Joc.Add(J0);

					Goc = new Matrix(1, 1);
					Goc[0, 0] = 1 / Alpha;

					X = X + L * Teta0;
				}
				// Другая ветка
				else
				{
					Poc = new Matrix(Joc.Count, 1);

					for (int i = 0; i < Joc.Count; i++)
					{
						CountOfGesseValue++;
						Poc[i, 0] = Gesse(X, Joc.ElementAt(i), J0);
					}

					//Рассчитываем Loc
					Loc = Goc * Poc * Lj0;

					var Doc = new Matrix(Joc.Count, Joc.Count);

					for (int i = 0; i < Joc.Count; i++)
					{
						for (int j = 0; j < Joc.Count; j++)
						{
							CountOfGesseValue++;
							Doc[i, j] = Gesse(X, Joc.ElementAt(i), Joc.ElementAt(j));
						}
					}

					//	Matrix Doc = Gesses(X).Submatrix(Joc.ToArray(), Joc.ToArray());





					//Matrix Doc = Gesses(X).Submatrix(Joc.ToArray(), Joc.ToArray());
					Matrix temp = Loc.Transpose() * Doc * Loc - Loc.Transpose() * Poc * Lj0 * 2;
					//Рассчитываем TetaAlpha]'=
					CountOfGesseValue++;
					double dj0j0 = Gesse(X, J0, J0);
					Alpha = (temp[0, 0] + dj0j0);

					if (Alpha < 0)
					{
						throw new Exception();
					}

					if (Math.Abs(Alpha) < 1e-8)
					{
						Alpha = 0;
					}


					for (int i = 0; i < Size; i++)
					{
						L[i, 0] = 0;
						//Teta[i, 0] = 0;
					}

					for (int i = 0; i < Joc.Count; i++)
					{
						L[Joc.ElementAt(i), 0] = Loc[i, 0];
					}
					L[J0, 0] = -Lj0;

					//Рассчитываем вектор Teta
					for (int i = 0; i < Size; i++)
					{
						if (L[i, 0] < 0)
						{
							Teta[i, 0] = (dmin[i] - X[i, 0]) / L[i, 0];
						}
						else if (L[i, 0] > 0)
						{
							Teta[i, 0] = (dmax[i] - X[i, 0]) / L[i, 0];
						}
						else
						{
							Teta[i, 0] = double.PositiveInfinity;
						}
					}

					TetaAlpha = Math.Abs(Deltaj0 / Alpha);
					CountOfFuncValue+=2;
					eta = (CalculateF(X) - CalculateF(X + L * TetaAlpha)) / 2;


					//Func<double, double> f = (teta) => CalculateF(X + L * teta);
					//var min = GoldCut(f, -2, 2, 1e-15);
					//TetaAlpha = min;


					//Выбираем Teta0
					TetaJ0 = Teta[J0, 0];

					Teta0 = TetaJ0;

					for (int i = 0; i < Size; i++)
					{
						if (Teta[i, 0] < Teta0)
						{
							Teta0 = Teta[i, 0];
							I0 = i;
						}
					}

					if (Teta0 > TetaAlpha)
					{
						Teta0 = TetaAlpha;
					}

					if (Teta0 == TetaJ0)
					{
						X = X + L * Teta0;
					}
					else if (Teta0 == TetaAlpha)
					{
						//Добавляем новый индекс в множество j0
						if (Joc.Contains(J0))
						{
							throw new Exception();
						}
						Joc.Add(J0);

						Matrix GocTemp = new Matrix(Goc.Rows + 1, Goc.Columns + 1);
						Matrix newGoc = new Matrix(Goc.Rows, Goc.Columns);
						newGoc = Goc + Loc * Loc.Transpose() * (1 / Alpha);

						Matrix underGoc = -Loc.Transpose() * (Lj0 / Alpha);
						GocTemp[GocTemp.Rows - 1, GocTemp.Columns - 1] = 1 / Alpha;
						Matrix rightGoc = -Loc * (Lj0 / Alpha);

						for (int i = 0; i < newGoc.Rows; i++)
							for (int j = 0; j < newGoc.Columns; j++)
							{
								GocTemp[i, j] = newGoc[i, j];
							}
						for (int i = 0; i < underGoc.Columns; i++)
						{
							GocTemp[GocTemp.Rows - 1, i] = underGoc[0, i];
						}
						for (int i = 0; i < rightGoc.Rows; i++)
						{
							GocTemp[i, GocTemp.Columns - 1] = rightGoc[i, 0];
						}
						GocTemp[GocTemp.Rows - 1, GocTemp.Columns - 1] = 1 / Alpha;

						Goc = GocTemp;
						X = X + L * Teta0;
					}
					else
					{

						var GocTemp = new Matrix(Goc.Rows - 1, Goc.Columns - 1);

						// Запоминаем индекс выводимого элемента
						var index = Joc.IndexOf(I0);

						//Запоминаем оставшиеся индексы


						var indexes = new List<int>();
						//foreach (var item in Joc)
						//{
						//	if (Joc.IndexOf(item) != index)
						//	{
						//		indexes.Add(Joc.IndexOf(item));
						//	}
						//}

						indexes.AddRange(Joc.Where(y => Joc.IndexOf(y) != y));

						//Выводим индекс из множества
						Joc.Remove(I0);

						//Формируем новую матрицу Goc с выведенным индексом
						for (int i = 0; i < Goc.Rows - 1; i++)
							for (int j = 0; j < Goc.Columns - 1; j++)
							{
								GocTemp[i, j] = Goc[indexes[i], indexes[j]];
							}

						Matrix GocHorizontal = new Matrix(1, Joc.Count);
						for (int i = 0; i < Joc.Count; i++)
						{
							GocHorizontal[0, i] = Goc[index, indexes[i]];
						}

						Matrix GocVertical = new Matrix(Joc.Count, 1);
						for (int i = 0; i < Joc.Count; i++)
						{
							GocVertical[i, 0] = Goc[indexes[i], index];
						}

						Goc = GocTemp - GocVertical * GocHorizontal * (1 / Goc[index, index]);

						//Пересчитываем вектор x
						X = X + L * Teta0;
						Console.WriteLine("F= {0}", CalculateF());
						//Console.ReadKey();
					}

				}
			}
		}
	}
}
